Week of 2025.07.30¶
Goal of this week:¶
- Create clean visulizations that show the current state of the infrasound work.
- Coordinate meeting with Fred + Mario + Eric to discuss the current state of things and figure out what needs to further be worked on.
- Get Fred's lightpath data to verify methods for calculating pressure at observer
- Put together slide show.
- Get access from to
/oak/stanford/groups/edunham
Things to change¶
- Plot mach numbers
- Color bar of atmosphere + plot heat
- Add atmosphere temperature plot.
- Reduce stem to maybe 10 m
- time series of flank of volcano.
Notes¶
- public wifi can block ssh access.
Sherlock notes¶
- Jupyter notebooks can be run via interactive sessions
- In theory there is a way to ssh onto sherlock via VS code like this. That appeals to me because it means I could run a simulation and visualize locally on vs code but have all the data live on sherlock.
- To run a job with dedicated resources do something like:
sh_dev -t 24which launches an interactive shell. - After running in some challenges with the sherlock time limit, I increased the time limit and changed the shell script to look like:
#!/usr/bin/bash #BATCH --job-name=quail #SBATCH --time=60:00:00 #SBATCH --ntasks=1 #SBATCH --cpus-per-task=16 #SBATCH --mem=32GB #SBATCH -p serc module load python/3.12.1 # User defineed variables export PATH=$PATH:/oak/stanford/groups/edunham/paxtonsc/quail_volcano/src cd $OAK/paxtonsc/volcano-sims/2025_07_21_eruption_sim/ quail conduit.py
This script I called run_quail.sbatch and I started the batch jobs with the call sbatch run_quail.sbatch and subsequently observed the state of the jobs with squeue --me.
Although as of this writing my job has not yet gotten to the time when it runs, I believe I can watch the progress with ssh sh02-01n49.
Questions¶
- Once you start a job on sherlock how the heck do you check its progress? Turns out a log is generated in the same folder as the bash script.
In [5]:
%load_ext autoreload
%autoreload 2
from helper_code.slip_imports import *
from helper_code.helper_functions import get_local_solver_from_index_func, get_quantities_at_conduit_exit, get_quantities_at_all_space
from helper_code.animate import animate_conduit_pressure, animate_melt_atmosphere_combination
from helper_code.analytical_solution import get_lumped_solution
import helper_code.infrasound as infrasound
from matplotlib.animation import FuncAnimation
import matplotlib.pyplot as plt
plt.rcParams['animation.ffmpeg_path'] = '/opt/homebrew/bin/ffmpeg'
folder_name = "simple_conduit"
ITERS = 50
D_ITERS = 1
END_TIME = 10
C0 = 320 # m/s at 5000m
The autoreload extension is already loaded. To reload it, use: %reload_ext autoreload
0.0 Velocity & Pressure at outlet of various runs¶
In [6]:
file_name = "tungurahua_no_atm_v0"
folder_name = "2025_07_21_eruption_sim"
solver_func = get_local_solver_from_index_func(folder_name, file_name)
t_tung_no_atm, p_tung_no_atm, slip_tung_no_atm, u_tung_no_atm, _ = get_quantities_at_conduit_exit(
solver_func, iterations=100, compute_temp=True, conduit_index=-1
)
file_name = "tungurahua_atm_1m"
folder_name = "2025_07_21_eruption_sim"
solver_func = get_local_solver_from_index_func(folder_name, file_name)
t_tung_atm_1m, p_tung_atm_1m, slip_tung_atm_1m, u_tung_atm_1m, _ = get_quantities_at_conduit_exit(
solver_func, iterations=100, compute_temp=True, conduit_index=-1
)
folder_name = "2025_07_23_eruption_1d_atmosphere_sim"
file_name = "tungurahua_atm_1m"
solver_func = get_local_solver_from_index_func(folder_name, file_name)
t_tung_atm_1m_1d, p_tung_atm_1m_1d, slip_tung_atm_1m_1d, u_tung_atm_1m_1d, _ = get_quantities_at_conduit_exit(
solver_func, iterations=100, compute_temp=True, conduit_index=-1
)
In [7]:
# Plot results
plt.figure(figsize=(10, 8))
plt.subplot(3, 1, 1)
plt.plot(t_tung_no_atm, slip_tung_no_atm, "r--", label="Eruption without atmosphere [1m mesh]")
plt.plot(t_tung_atm_1m, slip_tung_atm_1m, "b-", label="Eruption with atmosphere [1m mesh]")
plt.plot(t_tung_atm_1m_1d, slip_tung_atm_1m_1d, label="Eruption with atmosphere [1m mesh] [1D]")
plt.xlabel("Time [s]")
plt.ylabel("s [m]")
plt.grid(True)
plt.legend()
plt.subplot(3, 1, 2)
plt.plot(t_tung_no_atm, u_tung_no_atm, "r--", label="Eruption without atmosphere [1m mesh]")
plt.plot(t_tung_atm_1m, u_tung_atm_1m, "b-", label="Eruption with atmosphere [1m mesh]")
plt.plot(t_tung_atm_1m_1d, u_tung_atm_1m_1d, label="Eruption with atmosphere [1m mesh] [1D]")
plt.xlabel("Time [s]")
plt.ylabel(r"$\dot{s}$ [m/s]")
plt.grid(True)
plt.legend()
plt.subplot(3, 1, 3)
plt.plot(t_tung_no_atm, p_tung_no_atm, "r--", label="Eruption without atmosphere [1m mesh]")
plt.plot(t_tung_atm_1m, p_tung_atm_1m, "b-", label="Eruption with atmosphere [1m mesh]")
plt.plot(t_tung_atm_1m_1d, p_tung_atm_1m_1d, label="Eruption with atmosphere [1m mesh] [1D]")
plt.xlabel("Time [s]")
plt.ylabel("Pressure [Pa]")
plt.grid(True)
plt.legend()
plt.savefig(f"{BASE_PATH}/volcano_sims/notebooks/charts/fred_presentation/slip_pressure_slip_dot.png", dpi=300)
plt.tight_layout()
plt.show()
1.0 Eruption paired with 1D atmosphere¶
In [91]:
folder_name = "2025_07_23_eruption_1d_atmosphere_sim"
melt_file_name = "tungurahua_atm_1m"
atmosphere_file_name = "vertical_atmosphere_01_atm1"
melt_solver_func = get_local_solver_from_index_func(folder_name, melt_file_name)
atmosphere_solver_func = get_local_solver_from_index_func(folder_name, atmosphere_file_name)
iters = 100
ani = animate_melt_atmosphere_combination(
melt_solver_func,
atmosphere_solver_func,
iterations=100,
d_iterations=8,
y_min=-1000,
y_max=1000,
max_pressure=0.13,
min_pressure=0.07,
max_velocity=50,
min_velocity=0,
max_density=3e3,
min_density=0,
max_slip=150,
min_slip=0,
max_speed_of_sound=100
)
ani.save(
f"{BASE_PATH}/volcano_sims/notebooks/animations/fred_presentation/1d_atmosphere_combination.mp4",
writer="ffmpeg",
fps=10,
bitrate=1800,
)
HTML(ani.to_html5_video())
Out[91]:
In [58]:
file_name = "tungurahua_no_atm_v0"
folder_name = "2025_07_21_eruption_sim"
iters = 100
solver_func = get_local_solver_from_index_func(folder_name, file_name)
ani = animate_conduit_pressure(
solver_func,
iterations=iters,
d_iterations=1,
viscosity_index=1,
wall_friction_index=5,
max_velocity=20,
max_slip=60,
max_tau=1.5,
max_pressure=25,
max_speed_of_sound=2000,
max_water=1,
max_density=5e3,
max_fragmentation=5000,
max_mach=2.0,
max_viscosity=2,
)
ani.save(
f"{BASE_PATH}/volcano_sims/notebooks/animations/fred_presentation/conduit_flow_connected_2D_atmosphere.mp4",
writer="ffmpeg",
fps=10,
bitrate=1800,
)
HTML(ani.to_html5_video())
Out[58]:
2.2 Atmosphere¶
- Time elapsed for me to run this simulation locally: 10hr?
In [10]:
from helper_code.atmosphere_animate import atmosphere_2d_animate
folder_name = "2025_07_21_eruption_sim"
file_name = "test_infrasound_v1_atm1"
iterations = 80
ani = atmosphere_2d_animate(
iterations=iterations,
d_iter=4,
folder=folder_name,
file_name=file_name
)
ani.save(
f"{BASE_PATH}/volcano_sims/notebooks/animations/fred_presentation/2d_atmosphere.mp4",
writer="ffmpeg",
fps=10,
bitrate=1800,
)
HTML(ani.to_html5_video())
Out[10]:
2.1 Time Series Pressure Data¶
In [11]:
import matplotlib.tri as tri
from helper_code.infrasound import *
folder_name = "2025_07_21_eruption_sim"
file_name = "test_infrasound_v1_atm1"
solver_func = get_local_solver_from_index_func(folder_name, file_name)
iterations = 80
solver_0 = solver_func(0)
triangulation = tri.Triangulation(solver_0.mesh.node_coords[...,0],
solver_0.mesh.node_coords[...,1],
triangles=solver_0.mesh.elem_to_node_IDs)
trifinder = triangulation.get_trifinder()
p_20 = infrasound.get_p_series(20, 0, solver_func, trifinder, iterations=iterations, p0=1e5)
p_50 = infrasound.get_p_series(50, 0, solver_func, trifinder, iterations=iterations, p0=1e5)
p_80 = infrasound.get_p_series(80, 0, solver_func, trifinder, iterations=iterations, p0=1e5)
Element ID for point (20, 0): 5085 Element ID for point (50, 0): 6102 Element ID for point (80, 0): 6478
/Users/paxton/git/volcano_sims/notebooks/helper_code/infrasound.py:102: UserWarning: In TriCubicInterpolator initialization, PCG sparse solver did not converge after 1000 iterations. `geom` approximation is used instead of `min_E` U_target = np.array([tri.CubicTriInterpolator(
In [ ]:
fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(10, 6), sharex=True)
t_end = solver_func(iterations-1).time
t = np.linspace(0, t_end, len(p_20))
C0 = 320
axes[0].plot(t, p_20, label="pressure a r=20")
axes[1].plot(t, p_50, label="pressure at r=50")
axes[2].plot(t, p_80, label="pressure at r=80")
for i, ax in enumerate(axes):
ax.legend(loc='lower right')
ax.grid(True)
ax.set_ylim(-200, 200)
ax.set_xlim(0, t_end)
fig.suptitle("Pressure at various radi")
fig.supylabel("pressure (Pa)")
fig.supxlabel("time (s)")
plt.savefig(f"{BASE_PATH}/volcano_sims/notebooks/charts/fred_presentation/pressure_atm_time_series.png", dpi=300)
plt.tight_layout()
plt.show()
3.0 Light hill stress tensor¶
3.1 Accoustic pressure calculate as volume intergral¶
I am calculating the volume integral over the half-cylindrical region bounded by $r < 80 $ and $y > -80$ and $y<80$.
As we can see in the above graphics, the value of the double divergence of the lighthill stress tensor falls off outside those bounds.
The equation we solve below is:
$$ p'(x,t) = \frac{1}{4\pi} \int_v \frac{1}{|x - y|} \frac{\partial^2 T_{ij}}{\partial y_i \partial y_j} (y, t - \frac{|x-y|}{c_0}) dy^3 $$
To solve for pressure as a time series, we can just iterate over the integral for all time steps. We need only calculate the value of $\frac{\partial^2 T_{ij}}{\partial y_i \partial y_j}$ when $t - \frac{|x-y|}{c_0}$ is greater than 0 and less than the end of the simulation.
- The $4 \pi$ in the denominator is technically incorrect for our geometry.
My implementation can be found here. Validation for that implementation can be found in section 3.1 here.
In [97]:
import helper_code.lighthill as lighthill
import numpy as np
import matplotlib.tri as tri
iterations = 80
def solver_from_2D(dom, i):
solver_func = get_local_solver_from_index_func("2025_07_21_eruption_sim", f"test_infrasound_v1_atm{dom}")
return solver_func(i)
# Prep interpolation grid
solver0 = solver_from_2D(1, 0)
physics = solver0.physics
base_x = np.linspace(0, 125, 200)
base_y = np.linspace(-125, 125, 200)
mg_x, mg_y = np.meshgrid(base_x, base_y, indexing="xy")
solver0_list = [solver_from_2D(dom_idx, 0) for dom_idx in [1]]
# Compute workload partition
ind_partition = [np.where(tri.Triangulation(
solver.mesh.node_coords[...,0],
solver.mesh.node_coords[...,1],
triangles=solver.mesh.elem_to_node_IDs).get_trifinder()(mg_x.ravel(), mg_y.ravel()) != -1)[0]
for solver in solver0_list]
# List of file indices to read
file_index_list = np.arange(0,iterations,1)
In [98]:
# Allocate union (joining all times) U, in spatially-flattened shape
U_union = np.nan * np.empty((file_index_list.size, *mg_x.ravel().shape, 8+physics.NDIMS))
for time_idx, file_idx in enumerate(file_index_list):
# Load solvers for given time_idx
solver_list = [solver_from_2D(dom_idx, file_idx) for dom_idx in [1]]
for solver, _index_partition in zip(solver_list, ind_partition):
# Identify physical position (x, y) of points to interpolate at with shape (npoints, 2)
_phys_pos = np.stack(
(mg_x.ravel()[_index_partition],
mg_y.ravel()[_index_partition]),
axis=1)
# Identify element indices for all points to interpolate at
elt_indices = tri.Triangulation(
solver.mesh.node_coords[...,0],
solver.mesh.node_coords[...,1],
triangles=solver.mesh.elem_to_node_IDs).get_trifinder()(_phys_pos[:,0], _phys_pos[:,1])
# Identify element node coordinates
x_tri = solver.mesh.node_coords[solver.mesh.elem_to_node_IDs]
# Compute global physical-to-reference coordinate mapping
ref_mapping = lighthill.compute_ref_mapping(x_tri)
# Interpolate for each point using the correct element, writing to correct index in global U array
for (write_idx, x_point, ie) in zip(_index_partition, _phys_pos, elt_indices):
# For element containing point, compute reference coordinate of sampling point
ref_coords_loc = np.einsum("ij, j -> i",
ref_mapping[ie,...],
x_point - x_tri[ie,0,:])
# Evaluate basis at reference coordinate
#U_union[time_idx,write_idx,:] = (solver.state_coeffs[ie,0,:] * (1 - ref_coords_loc[0] - ref_coords_loc[1]))
# I believe this is the correct calculation for order 1
U_union[time_idx,write_idx,:] = (solver.state_coeffs[ie,0,:] * (1 - ref_coords_loc[0] - ref_coords_loc[1]) + solver.state_coeffs[ie,1,:] * ref_coords_loc[0] + solver.state_coeffs[ie,2,:] * ref_coords_loc[1])
if write_idx % 10 == 0:
print(f"Interpolated time index {time_idx}, file index {file_idx}, write index {write_idx} for element {ie}")